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We consider a single Josephson junction in the presence ol time varying gate charge, and examine 
the nonequilibrium work done by the charge control in the framework of fluctuation theorems. We 
obtain the probability distribution functions of the works performed by forward protocol and by 
its time reversed protocol, which from the Crooks relation gives the estimation of the free energy 
changes A_F = 0. The reliability of AF estimated from the Jarzynksi equality is crucially dependent 
on protocol parameters, while Bennett's acceptance ratio method confirms consistently AF = 0. 
The error of the Jarzynski estimator either grows or becomes saturated as the duration of the 
work protocol increases, which depends on the protocol rapidity determining the existence of the 
oscillatory motion of the phase difference across the junction. The average of the work also shows 
similar behaviors and its saturation value is given by the relative weight of the oscillatory trajectory 
with respect to running trajectories with constant acceleration. We discuss non-negativity of the 
work average and its relation to heat and entropy production associated with the circuit control. 

PACS numbers: 85.25.Cp, 05.70.Ln, 05.40.-a 



I. INTRODUCTION 

According to the first law of thermodynamics, the in- 
ternal energy change occurs due to the work (w) done 
to the system by external forces, and the heat input (Q) 
from the reservoir. Both w and Q are path dependent 
quantities but the internal energy change remains as a 
state function. Recently, Jarzynski found that the work 
alone satisfies a fluctuation theorem [l|. 
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Here the work is performed by applying a time dependent 
perturbation to a system. Let the control parameter of 
the perturbation, say \{t), vary in time t according to 
a prescribed protocol. The work defined by Jarzynski 
reads 



w 



d\{t) dH 



dt dX 
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where the integration interval corresponds to the time 
span of the work protocol. Initially the unperturbed sys- 
tem with A(0) is supposed to be in equilibrium state with 
a reservoir at the inverse temperature /3, and hence the 
initial phase space is populated according to the canoni- 
cal distribution e"^^^^^^'''' /Zi with the initial partition 
function Zi related to the initial free energy —/3Fi — 
InZi. Due to this randomness of the initial states, the 
work becomes a stochastic variable that can be character- 
ized by its probability distribution function p{w) which 
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can be acquired by many repetitions of the work pro- 
tocol. The left hand side of Eq. ([1]) is the average of 
the exponentiated negative work with respect to 
which, as the equality indicates, gives the free energy 
change AF = Ff — Fi where Ff denotes the free en- 
ergy of the final equilibrium state of the system described 
by the parameter A(r). The feasibility of this scheme 
have been demonstrated in various systems such as sin- 
gle molecules 0,11], classical oscillators and an elec- 
tronic circuit [7|. 

Also the protocol can be performed in bidirectional 
way. For the forward protocol choosing X{t), we call the 
work done during this process as the forward work, and 
write its probability distribution function as pf{w). For 
the backward protocol, the system is initially in thermal 
equilibrium with the perturbation at A(r), and the pa- 
rameter varies along the time-reversed path of \{t) and 
\{t) at the end of the backward protocol reaches A(0). 
We let the probability distribution function for this back- 
ward work to be pb{w). These two distributions, ppiw) 
and pb{w) are related through the Crooks relation @, 



Pf{w) = e-^'-'^^-'^^pBi-w). 



(3) 



This relation gives an identity pf{^F) = pb{—^F), 
indicating that the crossing point between pf(w) and 
Pb{—w) corresponds to AF, and hence allows to infer 
AF without the average process required in Eq. ([T]). 

In this work, we consider a system of a single Joseph- 
son junction and investigate the statistics of the work 
done by a gate charge control. Although the microscopic 
origin of involved mechanism in a Josephson junction is 
of quantum mechanical nature, certain features such as 
Josephson effects and Shapiro steps are explained purely 
on the basis of a classical equation of motion for the phase 
of the condensate wave function 0,113. We adopt this 
effective approach and show that the work performed by 
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a time varying gate charge satisfies the fluctuation the- 
orems ((H) and This proposes a possibility to verify 
the fluctuation theorems on a quantum mechanical sys- 
tem yet in the context of classical mechanics. Section HIl 
introduces our system and a work protocol to vary the 
gate charge in time. We consider bidirectional protocols 
which are related with each other via time reversal op- 
eration, and define the Jarzynski work for each direction 
of the work protocol. We also point out an equivalence 
to a current controlled junction. In Sec. Illli we obtain 
the probability distribution functions (pdfs) of the for- 
ward and the backward work, using Monte Carlo simula- 
tion along with the second order Runge-Kutta algorithm. 
The shapes of pdfs are markedly different depending on 
the protocol rapidity, and yet a symmetry relation be- 
tween the forward and the backward work exists to give 
AF — according to the Crooks relation. Sec. IIVI is de- 
voted to the free energy estimation by using the Jarzynski 
equality, the utility of which is limited by finite size of 
the work data. The error in the Jarzynski estimation of 
AF grows rapidly with the time elapse of the protocol, 
or becomes saturated. In comparison, the Bennett's ac- 
ceptance ratio method, which is an integrated form of 
the Crooks relation with a Fermi-function like weighting 
function introduced [ll|> is shown to give reliable esti- 
mations of AF irrespectively of the protocol parameters. 
Finally in Sec. |Vl we detail the properties of the work av- 
erage and its relation to heat and entropy generated by 
the circuit control, and summary will follow in Sec. IVII 

II. SYSTEM AND WORK 

We begin with the effective Hamiltonian for a single 
Josephson junction in the presence of a gate charge Qeit), 

H{t) - ^[q(t)-q,{t)f -EjCO^eit), (4) 

where represents the phase difference of condensed 
Cooper pairs between the two superconducting grains 
separated by an insulating barrier. The first term de- 
scribes a charging energy for excess charge q — qe and the 
second term denotes the Josephson coupling energy. We 
assume that the system is initially in equilibrium with 
a thermal reservoir at inverse temperature (5 so that the 
variables at initial time i = 0, g(0) and 0(0), are dis- 
tributed according to the canonical distribution 

Pe,(g(0),f?(0))=e-^^("V^. (5) 

with the initial partition function Zi — J d9 J dqe^^^'^^K 
For the classical picture, the equation of motion for 9{t) 
can be obtained from the Hamilton's equation: dtd{t) ~ 
{2e/h)dq'H and dtq = -{2e/h)deH with = d/dx, 
yielding 

9^0(0 = -2[sin0(<) -(- adtqe{t)l (6) 

where the ratio between the Josephson coupling en- 
ergy and the charging energy of a single Cooper pair 



Ec = (2e)2/(2C) enters as a = ^Ec/Ej. In Eq. dH) 
and from now on, the time t, the charge, the energy (and 
the work), and the inverse temperature are measured in 
units of Ti/ \/EcEj, 2e, Ej, and 1/Ej, respectively. In 
dimensionless form, the Hamiltonian Q is written as 

H{t) = a^[q{t) - q,{t)]^ - cos0(i). (7) 

Equation ^ equivalently describes the motion of a par- 
ticle moving in a tilted washboard potential. Depend- 
ing on the tilting slope, the particle can be either in the 
locked state staying around the potential minima or in 
the running state rolling down the hill [l2l | . In the present 
problem, the gate charge, qe{t), provides the tilting slope, 
and its time variance performs the nonequilibrium work 
on the Josephson junction. 

Let us in particular consider the simplest protocols for 
qe{t). For the forward process, we choose 

qe{t) = (7/a)t (8) 

and its time reversed path qf{t) = {'j/a){T — t) is taken 
for the backward process, where the both processes take 
place for a fixed time span t = {0,T). The parame- 
ter 7 in the protocol controls how rapid the gate charge 
changes in time. The Jarzynski work in our system can 
be obtained by taking A(i) ~ qe{t), which for the forward 
process reads 

w = -2aj [ dt[q{t) - q,{t)] - -j[9{T) - 9(0)], (9) 
Jo 

in units of the Josephson coupling energy Ej. For the 
second equality, we have used {l/2a)dt9{t) — q{t) — qe{t). 
For the backward process, denoting the angle path as 
9^{t), we write the work done by the time varying q^{t) 

WB^7K{T)-9^m, (10) 

in units of Ej. If the backward angle trajectory 9^{t) 
is the time reversal of the forward trajectory 9{t), that 
is, 9^{t) = 9{T — t), we expect a symmetry relation be- 
tween the forward work and the backward work to give 
Pf{w) — Pb{w)- Combining this with the Crooks rela- 
tion Eq. ([3]), we YiSive pf{w) = e~^^^ e^^pB{—w), which 
upon inserting w — Q gives the free energy change of this 
system AF — 0. This could be expected by that the 
classical partition function is invariant under the trans- 
lation of the momentum q{t) by qeit). Let us point out 
that the system described above can be equivalent to a 
current-biased situation with a special protocol for vary- 
ing the bias. In the presence of current bias, I{t) in units 
of {2e/h)^/EcEj, a Josephson junction can be described 
by the dimensionless Hamiltonian: 

H{t) = a^qitf - cos 9{t) + I{t)9{t), (11) 

which again leads to the equation of motion for 9 similar 
to Eq. dS]): 

d^9{t) = -2[sm9{t) + I{t)]. (12) 
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This implies a condition for an equivalence to the charge 
control: adtQeit) = I{t). For qe{t) — {'-f/a)t, one chooses 
a protocol I{t) = 7[e(t) - e(t - T)] with 8(2;) = 1 (0) 
for a: > (a; < 0) being the Heaviside step function, and 
then, Eq. (|T2t governing the time evolution of the phase 
angle 9(t) becomes identical to Eq. (|6|) with the initial 
distribution of the dynamic variables, 0(0) and q{Q), also 
determined by Eq. ([5]). Furthermore, the work in this 
case should be determined by taking X{t) — I{t), 

[ dte{t)dtl{t), (13) 
Jo 

which is again identical to the work by the gate control 
^ . In the foregoing discussion, we only refer to the case 
of gate charge control and results to be presented can 
directly be applicable to this current controlled junction. 



III. PROBABILITY DISTRIBUTION 
FUNCTION OF THE WORK AND THE CROOKS 
RELATION 

In order to obtain the pdfs of the work defined by 
Eqs. (O and ^TU\i for the forward and the backward pro- 
tocols, we first generate the initial dynamic variables, 
6{0) and q{0), according to the canonical distribution 
Eq. ([5]). Provided these initial values, the equation of 
motion © is solved through the second order Runge- 
Kutta method. We acquire 10^ work values each for the 
forward and the backward protocol, and present their 
distributions in Fig. [TJ We find that the resulting distri- 
bution of the forward work, pp{w), and the distribution 
of the negative of the backward work, pb{—w) are mir- 
ror symmetric about w = 0, illustrating the symmetry 
property, pf{w) = pb(w). As anticipated, they cross 
each other at = 0, which from the Crooks relation 
^ indicates AF — 0. In other words, the perturba- 
tion leaves the system free energy unaltered. This result 
should not depend on the protocol rapidity 7 in Eq. (|8]) 
which only changes the shape of the distributions. At 
7 = 2 [see Fig. HJa)], the pdf of the work appears to be 
of the Gaussian form. For large 7, the nonlinear force 
term sai9{t) in Eq. ([5]) can be neglected in comparison 
with the second term given by 7 for the protocol (|S]). For 
this case, the phase trajectory is approximately given by 
e{t) = 6'(0) -I- VQt + -it^ with vq = dtO{t)\t=Q. This run- 
ning state yields the forward work defined by Eq. ©, 
w = —^{vqT + 7T^). The average and the second cum- 
mulant of this work are obtained as 

{w) = {-lTf = Wr (14) 
cj^ ^ {w') - = 2(7r)V/3. 

In Fig. [Ija) , we present the Gaussian distribution func- 
tion given by P(u') = exp[-{w~Wr)^ /{2a'^)]/V2TTa'^ (the 
curve with shaded area), which coincides well with the 
numerically obtained pf{w). 
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FIG. 1: Probability distribution functions (pdfs) of w and 
—WB (in units of Ej), defined in Eqs. ([9]) and (fTO|) . respec- 
tively, where we choose the temperature for the initial thermal 
equilibrium state of the system, P — 0.1 (in units of 1/Ej), 
and the protocol duration time, 'yT = 3\/T0. The pdf of the 
forward work, ppiw), is mirror symmetric with the pdf of 
the negative of the backward work, pb{—w), and their cross- 
ing point at u) = gives AF — from the Crooks relation. 
Panel (a) presents the case for j — 2, where the work val- 
ues are almost normally distributed, as indicated in the co- 
incidence with the Gaussian distribution function (the curve 
with shade) with average and the variance given by Eq. (|14|l . 
Panel (b) displays the work distributions for 7 — 0.2 with a 
sharp peak pronounced around w = due to the oscillatory 
motion of the phase angle. Note that the side peaks are al- 
most identical to the distribution for 7 = 2 apart from the 
skewness (here the shaded region is identical to that presented 
in the panel (a)), indicating its origin from the running state 
agitated at finite temperatures. 

On the other hand, when 7 is small, there exits a 
locked state around a force-free point, 0^ ~ ~ sin~^ 7. 
Linearizing Eq. ^ around 6c, we obtain an oscillatory 
solution, 0{t) ~ 9c — Acos{ilt + (j>), where the amplitude 
A and the initial phase ((> are determined by initial con- 
ditions, 6'(0) - 6'c = A cos and dtO{t)\t=o = 2ag(0) = 
—AQsincj}. Here the oscillation frequency, fl is given 
hy il^ = 2| cos 6'c I . This oscillatory trajectory yields the 
work average, 

(w) = jec[cos{nT) -l]=wi , (15) 

which oscillates in the protocol duration time T with the 
amplitude much smaller than the average from the run- 
ning state [see Eq. for T ^ 1. As displayed in 
Fig. [TJb) for 7 = 0.2, this oscillation results in the sharp 
peak around w « in the pdf. Here, the broad side 
peak has large overlap with the Gaussian distribution of 
the running state, indicating that it originates from the 
running motion for initial high energy state at a finite 
temperature. 
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IV. FREE ENERGY ESTIMATION FROM THE 
JARZYNSKI EQUALITY AND THE BENNETT 
METHOD 

We have shown that the probabihty distribution func- 
tions Pf{w) and p_b(— uj) obtained from the bidirectional 
protocols cross with each other at w = /S.F — 0, illus- 
trating the feasibility of the Crooks relation in deter- 
mining the free energy change caused by the time vary- 
ing gate charge in a single Josephson junction. This 
should be confirmed also from the Jarzynski identity. 
Within the limitation on the number of work measure- 
ments {N = 10^ in this study), the estimation of free 
energy difference IS.F can be made from 



/3AFf = - 
P^Fb = In 



In 



1 ^ 



-/3tOi 



(16) 



1 ^ 
N ^ 



where wi and Wi^B denote the ith realization of the work 
in the forward and the backward protocol, respectively. 
Since the backward work is symmetrically distributed 
with respect to the pdf of the forward work, as exem- 
plified in Fig. m we have I^Fb ~ — AF, and finite AF 
indicates the bias of the estimation due to the finiteness 
of sample size TV. 

It is often the case that the Jarzynski equality fails to 
give a reasonable estimation of AF because the exponen- 
tial work average crucially depends on the sa mpling of the 
tail part the work distribution with w ^ AF The 
finite sampling error in the large sampling limit was in- 
vestigated in several studies (iSl - llSt . Probing the region 
w ^ AF is also important in extracting AF from the 
Crooks relation for the crossing criterion between priw) 
and pb{~w). In our case, the exact value of free en- 
ergy change is AF = 0, and therefore, the reliability of 
the Jarzynski identity and the Crooks relation depends 
on how low the probability to observe the work w ^ 0. 
For the Gaussian distribution displayed in Fig. [T] (see the 
upper panel), upon increasing 7T, the center of the dis- 
tribution given by the average (jl4p moves toward large 
positive values, and p{w ^ 0) becomes extremely small. 
In this case, the Jarzynski estimation can possess a severe 
error, and also the Crooks relation likely fails to give AF 
due to the absence of the overlap region between pp(w) 
and pb{—w). 

It is well worth introducing another method which is 
less restrictive than those, called the Bennett's acce 
tance ratio method based on the following equation 



cep - 



dw 



ppiw) 



dw 



Pb{~w) 

I ^ g/3(A_F-«,) ' 



(17) 



which can be obtained by multiplying the both sides 
of Crook's relation Q by a Fermi function 1/(1 -|- 
g/3(u;-A_F)'j_ 'pjjjg equation was originally implemented in 
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FIG. 2: (a) The bias of the free energy change calculated by 
using Eq. (|16p and the Bennett's acceptance ratio method as a 
function of the protocol parameter, 'yT, at the inverse temper- 
ature j3 — 0.1 (in units of 1/Ej). The Bennett method (line 
for 7 = 0.2 and triangles for 7 = 2) persistently yields 
AF ~ 0. but the Jarzynski estimation gives incorrect val- 
ues of AF for 7T > 10. The error for 7 = 2 (filled circles) 
rapidly grows as 'yT increases, whereas for 7 — 0.2 (open 
squares) the estimated AF becomes saturated. Panel (b) 
details this saturation behavior, which also shows the tem- 
perature dependence; the saturated error is smaller at lower 
temperature (/3 = 0.5). 



order to estimate partition function ratios by means of 
Monte Carlo sampling, where the Fermi-function weight 
in Eq. (I17p was derived from the requirement of mini- 
mal variance of the partition function estimator in the 
large sample limit. Recently, it was shown by Shirts et 
al. fl^ that a maximum likelihood estimate of the free 
energy change also yields that Bennett's acceptance ratio 
method. In various studies, the estimation of AF from 
Eq. (|17p was demonstrated to outperform the Jarzynski 
method and Crooks' crossing criterion [20l - |23[ . For an 
equal size N of forward and backward samplings, Eq. p7)) 
is written as 



1 



N 
i—1 



N 



1 



(18) 



which yields the free energy estimation AF from the Ben- 
nett's acceptance ratio method. In practice, especially 
when 7 is large we observe that the values of the left- and 
the right-hand sides of Eq. ([T8|) become extremely small 
below the numerical accuracy of computer in a broad 
range of AF. In this case, a precise estimation of AF 
from Eq. (|18p is not plausible, and we instead approxi- 



5 



mately estimate AF as {AFf-AFb)/2 \2^, where AFp 
and AFb are as defined in Eq. (|16p . 

In Fig. [21 we present results for tlie free energy es- 
timations from the Jarzynski relation and the Bennett 
method. As shown in the upper panel of Fig. [21 es- 
timations based on the Bennett method are very close 
to the true value {AF = 0) almost insensitively of the 
protocol parameter, jT. Meanwhile, the Jarzynski es- 
timations have finite errors. For 7 = 2 [see the filled 
circles in Fig. [Ha)], in particular, the error becomes 
more significant as increases. In fact, in order to 
obtain a relatively reliable estimation of AF from the 
Jarzynski equality, the number of acquired work values 
should be larger than Nc ~ e^'' with the hysteresis, 
h = {{w) + {wb))/^ From this one can estimate 

the onset point of the error by solving N = Nc, 



/3(7T)2=lniV, 



(19) 



where we used h k, Wr — {iTY. For the sample size, 
N = 10^, and the inverse temperature /? = 0.1, we find 
{^T)c ~ 12 at which indeed, the Jarzynski error for 7 = 2 
in Fig. begins to rise rapidly. Meanwhile, the error 
for 7 = 0.2 [see the open squares in Fig. [H^a)] is much 
smaller than for 7 = 2, and shows a saturation. As re- 
marked, when 7 < 1, there exists an localized motion 
of the phase variable near the potential minima which 
yields the work w fa 0, as manifested by the sharp peak 
at the origin of the pdf in Fig. [2jb) . These work values 
make the dominant contribution to Eq. (jl6p where the 
contributions from positively large w populated in the 
side broad peak for the running motions are exponen- 
tially small. These work values pinned near the origin are 
responsible for the smaller bias than the case of 7 > 2, 
and also for the saturation. Roughly, we take the work 
values from the localized oscillatory motions to be w = 
and let Ni to be the number of such realizations of the 
work values. Then we have 



13 AF « - lii{Ni/N), 



(20) 



which tells that the saturated value of AF is the frac- 
tion of the localized trajectories relative to the running 
motions. This fraction depends on the temperature, 
as shown in Fig. [Ifb), where at the lower temperature 
/3 = 0.5, the saturated value of AF becomes smaller, 
well corresponding to an expectation that the running 
state must be more suppressed at low temperature. The 
relative fraction between the localized and the running 
motion comes into play also in the work average, which 
will be discussed in the next section. 



V. WORK AVERAGE AND ENTROPY 
PRODUCTION 

Another important consequence of the fiuctuation the- 
orems is the second law of thermodynamics. The Jarzyn- 
ski equality, {e~^^) = e~^^^ together with the Jensen's 
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FIG. 3: (a) The work average is depicted as a function of the 
total time elapse T of the work protocol at various temper- 
atures, where the ordinate is chosen to be the work average 
divided by (7^)^ = Wr in Eq. (|14p corresponding to the work 
average by the running state. The values for 7 = 2 (filled sym- 
bols) are very close to unity, showing the dominance of the 
running state, at any temperatures presented here (tempera- 
tures are denoted in the lines of corresponding empty symbols 
for 7 = 0.2). For 7 = 0.2, the work average appears to be 
saturated into a certain fraction of Wr which becomes reduced 
by lowering temperature (i.e., increasing /3). Panel (b) shows 
the oscillation of the work average for 7 = 0.2 at extremely 
low temperature /3 = 10, which reaches the minimum value 
when the protocol duration T is close to an integer multiple 
of the phase oscillation period, 2% /Q,. 



inequality, e~''<'"> < (e"^'") lead to (w) > AF. Let 
us remind of relations, AF = AU — AS/{kBP) and 
{w) — AU = (Q). Here Q is the heat absorbed by the 
reservoir during the equilibration process; the system at 
the end of the work protocol is in nonequilibrium state 
and can reach the final equilibrium sufficiently long af- 
ter it is brought into a thermal contact with a reservoir. 
Then, the inequality, (w) — AF > can be translated 
into AS" + ASr > if the reservoir absorbs the heat Q 
reversibly from the system so that the entropy change of 
the reservoir is given by ASr = kB(3{Q) [H, [2^. Hence, 
the total entropy production in the system plus reser- 
voir is always positive: AStot — AS + ASr > 0. For 
AF = AU = as in our case studied here, the work 
average, identical to the average of the heat, determines 
the total entropy production: 



(w) - (Q) - AStot/ikBl3) > 



(21) 



Note here that the positivity of (w) for the second law of 
thermodynamics also indicates the direction of the heat 
from the system into the reservoir. 
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Having this in mind, let us examine the work average. 
As noted in Sec. IIV[ the work average crucially depends 
on 7. For 7 > 1, 9(t) increases in time, so does the 
work average. If nonlinear effect due to the sinusoidal 
force can be neglected for 7 3> 1, Wr for the running 
state given in Eq. (jl4[) . can be a good approximation 
of the work average. In Fig. |31[a), we present the work 
average for 7 = 2 as a function of the time elapse of 
the work protocol. In order to see more clearly the con- 
vergence behavior of (w) into = (7T)^, we choose 
the ordinate to be {w)/{'-fT)^. On the other hand, for 
7 < 1, in addition to the running state, the phase mo- 
tion can also be localized around Oc = — syoT^ 7 and 
yields the work average oscillating in time, as discussed 
in Sec. IIIII In this case, the work average can roughly 
be (w) sa {Nr/N)wr + {Ne/N)we with Nr = N - Ng. 
In particular for 7T 3> 1, since Wr ^ W£, we have that 
(w) « {Nr/N)wr- Therefore, the saturated value of the 
work averages shown in Fig.[3Jb) for 7 — 0.2 gives Nr/N. 
At (3 — 0.1, we have Nr/N 0.8, which is consistent with 
Ne w 0.2 estimated from the saturation of AF « 1.5 [see 
Fig. mb)] given in Eq. (120)) . This fraction of the running 
state decreases as lowering the temperature. At /3 = 0.5, 
we have Nr/N 0.52, well corresponding to AF = 0.74 
[see Fig. ^h)] giving Ne/N « 0.47. At the extremely 
low temperature, the contribution from the running state 
becomes vanishingly small and therefore, the oscillating 
behavior becomes clearly visible, as displayed in Fig.[3ljb) 
for /3 = 10. Recalling here the relation Eq. (PT|) stating 
that the work average is identical to the amount of the 
heat transferred from the system into the reservoir, and 
the entropy production, one finds that the heat and en- 
tropy generation associated with the circuit control cru- 
cially depends on 7 which determines the protocol ra- 
pidity by dtqe{t) — 7/a. For a fixed value of a and the 
terminal value of the gate charge qf = qe{T) = ^T/a, 
the heat generation for a fast protocol (7 > 1) is given 
by (Q) ^Wr = a'^qj = {Ec/Ej)qj (in units of Ej), cor- 
responding to the energy required to charge the junction 
of the capacitance C with the terminal gating charge g/. 
For a slow protocol giving 7 < 1, the heat generation in 
the large 7T limit can be expressed as (Q) ~ 
with an effective gate charge qf = \J Nr/Nqf modified 



by the temperature dependent factor Nr/N ^ which can 
be significantly suppressed by lowering temperature. 



VI. SUMMARY 

In summary, we considered a single Josephson junction 
and examine the nonequilibrium work done by a time 
varying gate charge. We found that the details of the 
probability distribution function of the work depends on 
the rapidity of the protocol. For a rapid protocol (7 > 1), 
the running state is dominant to yield the work values 
which are normally distributed. A slow protocol (7 < 1) 
allows the locked state which manifests itself by a pdf 
with sharp peak near the origin, which at finite temper- 
ature becomes bimodal with a broad side peak from the 
running state. A symmetry relation between the forward 
and the backward distribution, Pf{w) — pb{w), exist to 
yield AF = according to the Crooks relation. The 
Jarzynski estimation of AF can be biased for the finite 
sample size of the work data. The error grows rapidly 
with the time elapse if 7 > 1, and for a slow protocol 
7 < 1 the allowed oscillatory motion saturates the error 
into a certain value which is determined by the relative 
population of the locked state with respect to the running 
state. The Bennett's acceptance ratio method gives re- 
liable estimation, AF 0, irrespectively of the protocol 
parameters. We also discuss the behaviors of the work 
average reflecting the two characteristic motions, in re- 
lation to the heat and the entropy production associated 
with the charge control. 
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